In the field of medical imaging, various systems have been developed for generating medical images of various anatomical structures of individuals for the purpose of screening and evaluating medical conditions. These imaging systems include, for example, computed tomography (CT) imaging systems, magnetic resonance imaging (MRI) systems, X-ray systems, ultrasound systems, positron emission tomography (PET) systems, and single photon emission computed tomography (SPECT) systems, etc. Each imaging modality may provide unique advantages over other modalities for diagnosing and monitoring diseases, medical conditions or anatomical abnormalities. Three-dimensional modeling of image data has led to significant advances in both research and clinical capabilities.
A digital image consists of a two-dimensional (2D) array of data elements, pixels, that represent color or light intensity. Volumetric data are three-dimensional (3D) in structure, with a voxel as the smallest element. For example, computed tomography (CT) imaging systems can be used to acquire a volumetric dataset of cross-sectional images or two-dimensional (2D) “slices” that cover an anatomical part, allowing a 3D visualization of the part to be constructed.
Characterization of complex shapes embedded within volumetric data is an important step in a wide range of applications. Standard approaches to this problem employ surface based methods that require inefficient, time consuming, and error prone steps of surface segmentation and inflation to satisfy the uniqueness or stability of subsequent surface fitting algorithms. For example, the identification of structural changes in the brain on magnetic resonance imaging (MRI) scans is increasingly important in the study of neurological and psychiatric diseases. MRI can be used to identify and exclude treatable causes of cognitive impairment and it has also become important in the differential diagnosis of disease, in tracking disease progression, and for research purposes. Pathological changes in the brain resulting in cell loss manifest as loss of brain tissue, or atrophy, which can be detected by structural MRI. Characteristic patterns of atrophy are associated with specific neurodegenerative diseases. Traditional techniques of analyzing atrophy on MRI include visual assessment by experienced radiologists and manual measurements of structures of interest. However, automated techniques have been developed which allow the assessment of atrophy across large groups of subjects without the need for time-consuming manual measurements or subjective visual assessments.
Voxel-based morphometry (VBM) is one such automated technique that has grown in popularity since its introduction, largely because of the fact that it is relatively easy to use and has provided biologically plausible results. It uses statistics to identify differences in brain anatomy between groups of subjects, which in turn can be used to infer the presence of atrophy or, less commonly, tissue expansion in subjects with disease. The technique typically uses T1-weighted volumetric MRI scans and essentially performs statistical tests across all voxels in the image to identify volume differences between groups. For example, to identify differences in patterns of regional anatomy between groups of subjects, a series of t tests can be performed at every voxel in the image. Regression analyses can also be performed across voxels to assess neuroanatomical correlates of cognitive or behavioral deficits. The technique has been applied to a number of different disorders, including neurodegenerative diseases, e.g., Alzheimer's disease and dementia due to brain atrophy, movement disorders, epilepsy, multiple sclerosis, and schizophrenia, contributing to the understanding of how the brain changes in these disorders and how brain changes relate to characteristic clinical features. Although results from VBM studies are generally difficult to validate, studies have compared results of VBM analyses to manual and visual measurements of particular structures and have shown relatively good correspondence between the techniques, providing some confidence in the biological validity of VBM.
Standard approaches to volumetric data analysis employ surface based methods that require an initial segmentation of a surface and often a subsequent inflation of this surface to satisfy the uniqueness or stability of subsequent surface fitting algorithms. These methods are inefficient and time consuming because of the need for segmentation prior to fitting and the computationally intensive inflation process, the latter of which can be a significant source of errors due to creation of topological defects.
Continuing progress in the development of advanced imaging methods such as MRI and CT have facilitated the acquisition of high resolution, high contrast volumetric data that provides an approach for non-invasive yet highly informative assessment of brain morphology. The ability to utilize these valuable diagnostic and analytical tools in a cost effective manner is impeded by the massive quantity of data involved in modern volume images. This is particularly true of MRI data, which has a wide range of contrast mechanisms by which it can produce very high contrast of different types between complex soft tissues.
In concert with these advancements in imaging technologies, advances in computational methods, particularly in volume graphics and computer vision has resulted in tremendous increase in computational methods for morphology characterization and segmentation for comparative morphometry for both basic neuroscience studies on comparative brain anatomy and clinical studies of disease characterization and progression in humans, and for a broad range of studies in comparative biology.
In comparative biology, geometric morphometrics has emerged as an important tool for analysis, becoming commonly used to quantify morphology, wherein landmark points are identified in photographic (2D) images and then are fit to a warped mesh that provides a common coordinate system in which different specimens can be compared (Zelditch et al., 2004). These methods allow users to define key points of known morphological interest and statistically compare morphologies based on these points. However, the current predominant methods are based on 2D digital images or on 3D surfaces and are not readily applied to volumetric 3D data, such as those acquired by MRI or CT.
Recent advances in segmentation techniques were mostly originated from fuzzy logic and supervised and non-supervised clustering (Barra and Boire, 2000; Lin et al., 2012) both in 2D (Barra and Boire, 2001; Cocuzzo et al., 2011; Pedoia and Binaghi, 2012; Razlighi et al., 2012; Suri, 2001; Zavaljevski et al., 2000) and 3D (He et al., 2011; Kiebel et al., 2000; Klauschen et al., 2009; Popuri et al., 2012; Wels et al., 2011). Unfortunately, in spite of all advances none of these methods are able to provide truly robust and automated segmentation.
The most straightforward approach to segmentation is thresholding, which involves finding an intensity value, “the threshold”, which distinguishes features of interest. This method is most frequently used to create a binary segmentation of an image, but it is also possible to distinguish three or more intensity classes using multithresholding (Zavaljevski et al., 2000). Thresholding works particularly well with imaging modalities such as CT data where images are often essentially binary between bone (bright) and soft tissue (very dark) and segmentation can be practically automated. Automated methods for MRI data, however, are exceedingly difficult because of adjoining regions with similar values (i.e. low contrast), partial voluming (multiple tissue types within a single voxel), image noise, and intensity inhomogeneities, all of which are common to MR images (Atkins and Mackiewich, 2000; Pham et al., 2000).
Region growing methods extract connected regions in images based on criteria that can include both intensity and edges. These methods are susceptible to noise, which can create artificial divisions between connected regions, and partial volume effects, which can merge disconnected regions. These effects can be reduced by limiting growth to topology-preserving deformations (Mangin et al., 1995), but user input is still required to select seed regions. Clustering algorithms alternate between segmenting the image and characterizing the properties of each segmented class, iterating until a stopping criterion is reached (Barra and Boire, 2000, 2002; Liang et al., 1994; Pachai et al., 2001; Popuri et al., 2012). Clustering is generally susceptible to both noise and image inhomogeneities, though robustness to intensity inhomogeneities has been demonstrated (Pham and Prince, 1999). Given a Bayesian prior model, Markov random field models can be incorporated in clustering methods to minimize susceptibility to noise (Li, 1994).
Atlas-guided approaches provide an option that may be feasible (Klein et al., 2009). In such methods, a linear or non-linear transformation is found mapping the pre-segmented atlas image to the target image. This changes the tissue classification problem to a registration or deformation problem. However, to effectively use atlas-guided methods, very large and detailed databases, or atlases of reference objects, are needed. This puts the onus of the quantitation on an accurate and reproducible method for atlas creation.
One important and rather successful direction in brain quantifying and characterization has emerged from analyses of parameterization of surfaces for 3D shape description using spherical harmonics (SPHARM) representation (Brechbiihler et al., 1995; Kazhdan et al., 2003). SPHARM is a technique for parameterizing anatomical boundaries using the spherical harmonic basis. The surface coordinates are expressed as a linear combination of basis functions. SPHARM can be used to construct the Fourier series expansion of a functional measurement. Shape signatures can be created using the SPHARM decomposition at several concentric spheres or at a single surface that represents a highly convoluted geometry of the cerebral cortex. Two steps are involved in converting the object surface to its SPHARM shape description: surface parameterization and expansion into a complete set of SPHARM basis functions. In the SPHARM method, any function ƒ is assumed to be defined on a sphere, ƒ(θ, ϕ), and decomposed as the sum of its spherical harmonics:
                              f          ⁡                      (                          θ              ,              ϕ                        )                          =                              ∑                          l              =              0                        ∞                    ⁢                                          ⁢                                    ∑                              m                =                                  -                  l                                            l                        ⁢                                                  ⁢                                          f                lm                            ⁢                                                Y                  l                  m                                ⁡                                  (                                      θ                    ,                    ϕ                                    )                                                                                        (        1        )            where Ylm is the spherical harmonic of degree l and order m, with low values of l corresponding to lower frequency information. The Fourier coefficients ƒlm, are used to reconstruct the object surface, with a greater number of coefficients providing a more detailed construction. Since L2-norms of spherical functions are not affected by rotations, a rotationally-invariant shape signature may be given asSH(ƒ)={∥ƒ0(θ,ϕ)∥,∥ƒ1(θ,ϕ)∥, . . . }, whereƒl(θ, ϕ)=Σm=−ll ƒlmYlm(θ, ϕ) are the frequency components of ƒ. An alternate signature can be calculated more quickly and directly from the coefficients, defining SH(ƒ)={A0, A1, . . . }, where the Al are L2-norms of all the coefficients ƒlm at each l: Al=Σm=−ll|ƒlm|2. The spherical harmonics Ylm are continuous functions, but for computational applications, ƒ is only sampled at NΩ discrete angles. To create a shape signature for a 3D object, the shape is sampled at NΩ angles and Nγ radii, SH is calculated at each radius up to 1=Lmax, and the result is represented as a 2D grid of size Lmax×Nγ.
The SPHARM application described by Kazhdan et al. (2003) provided more general shape classification using clean data (e.g., a set of 1,890 household objects), but in noisy MRI data, the SPHARM deals with noise automatically, since the noise does not appear in the lower frequencies that dominate shape descriptions. Many internal structures remain visible in data reconstructed from the signatures, while the signatures themselves require significantly less storage space than the original data. This general method was improved further by appropriate filtering (i.e. using exponentially weighted Fourier or spherical harmonics series, Chung et al. 2008a, 2007, 2008b). The weighting reduces a substantial amount of the so called ringing (or Gibbs) phenomenon and aliasing (or Moire) patterns (Gelb, 1997)), both appearing because of relatively slow convergence of Fourier series when used for representing discontinuous or rapidly changing measurements.
Overall, these modifications of the SPHARM method with filtering or exponential weighting (Chung et al., 2008a, 2007, 2008b) allowed successful parameterization of the cortical surface including characterization of the local difference in gray matter concentration. Nevertheless, techniques based on the SPHARM morphometry method—tensor-based morphometry—uses the cortical surface already segmented out of noisy MRI data and quantifies the amount of gray matter only in a narrow layer along the tangential direction of the surface via the concept of a local area element. Hence, the analysis cannot be directly used for volumetric MRI data.
An extension of spherical harmonics decomposition that naturally allows incorporation of complete 3D volume data is known in various areas of physics, e.g., in quantum physics for description of waveform solutions of the Schrödinger equation (Gersten, 1971), in atomic and nuclear physics for approximation of Coulomb scattering function (Barnett, 1996, 1981)), and in astrophysics for analyses of anisotropies of microwave background, as well as for quantum gravity (Abbott and Schaefer, 1986; Binney and Quinn, 1991; Leistedt et al., 2012). However, such approaches have not previously been utilized for characterization and modeling of volumetric data.